Stationary states in single-well potentials under symmetric Levy noises 
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We discuss the existence of stationary states for subharmonic potentials V{x) oc |a::|'^, c < 2, under 
action of symmetric a-stable noises. We show analytically that the necessary condition for the 
existence of the steady state is c > 2 — a. These states are characterized by heavy-tailed probability 
density functions which decay as P{x) oc x~^'^'^'^~^' for |x-| — >■ oo, i.e. stationary states posses a 
heavier tail than the corresponding a-stable law. Monte Carlo simulations confirm the existence of 
such stationary states and the form of the tails of corresponding probability densities. 
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I. INTRODUCTION 

The dynamics of many physical systems can be reduced 
to a standard model of a "particle" (relevant coordinate 
x) under action of a deterministic force f{x) and the 
noise force including the action of all neglected degrees 
of freedom. The mathematical tool for such a description 
is given by a Langevin equation, which in the overdamped 
limit typically takes the form 



x{t) = j{x)+at). 



(1) 



In Eq. ^ f{x) stands for the deterministic force, while 
C(t) is the stochastic force — noise — describing inter- 
actions of a test particle with its complex surrounding. 
Commonly, it is assumed that the stochastic force is of 
the white type, for example representing the particle's 
collisions with molecules of the bath. Following General- 
ized Central Limit Theorem, a large number of indepen- 
dent collisions will lead to white noise of the stable type, 
i.e. either to Gaussian noise or to more general Levy 
noise. The presence of fluctuations distributed according 
to Levy laws have been observed in various situations 
in physics, chemistry or biology [l|, Q, paleoclimatology 
Q or economics [j|. The situations pertinent to heavy- 
tailed distributions appear in context of different models 
and are analyzed in an increasing number of studies 
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The action of white Levy noise leads to increments 
of the stochastic process which arc distributed accord- 
ing to a-stable Levy type distributions. Symmetric Levy 
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distributions Pa{x] a) are characterized by their Fourier- 
transforms (characteristic functions of the distributions) 
,/,(fc) = /^ e''=^p„(a;; a)dx being [H-il 



(^(/fc)=exp[-fT"|fc|"]. 



(2) 



The parameter a (where a € (0,2]) is the stability index 
of the distribution describing, for a < 2, the asymptotic 
decay of its tails 



p^ix) (X |xr(i+"). 



(3) 



Finally, a is the scale parameter which controls the over- 
all distribution width. The Gaussian distribution corre- 
sponds to a special case of a Levy stable law with a = 2 
when ^/2a is interpreted as the standard deviation of the 
distribution. In more general cases, for a < 2, the vari- 
ance of a-stable densities diverges. For a < 1, the mean 
value of a-stable densities also does not exist. 

The present work addresses properties of Levy flights 
in external potentials. Theoretical descriptions of such 
systems is based on the Langevin equation and/or 
Fokker-Planck equation, which is of the fractional order 



The research performed here extends earlier studies 
21| where analysis of Levy flights in harmonic and 
superharmonic potentials has been presented. The dis- 
cussion conducted there did not account for the problem 
of existence of stationary states in potential less steep 
than parabolic, which is the main topic of the present 
work. This issue is addressed here by the use of analyti- 
cal arguments and Monte Carlo simulations. 

The model under discussion is presented in Section HIl 
Section|lIl]discusses obtained results. The paper is closed 
with concluding remarks (Section IIV[) . 







FIG. 1. Exemplary subharmonic potentials: V{x) = \x\ and 
V{x) = |2:|^'^ used for examination of the problem of existence 
of stationary states. 



II. MODEL 

An overdampcd Brownian-Levy particle moves in an 
external potential V{x), therefore Eq. ([T|) takes the form 



x{t)^-v'{x) + at), 



(4) 



where V{x) represents a subharmonic single well poten- 
tial V{x) = Ix]" with < c < 2, see Fig. [Tl and C(t) de- 
notes a Levy stable white noise process [3, llflj UM, ISa ^M ■ 
The position of a random walker can be calculated by 
means of the stochastic integration of Eq. (j4]) [25|, |30] 



x{t)=x{0)- / V'{x{s))ds+ / C{s)ds 



x(0) 



V'{x{s))ds + La{t). 



(5) 



The integral ffj Cis)ds = 
process La{t) [3, [ifl, \ll 



Lnit) defines an a-stable Levy 
I l28l . |31| which is driven by a 
Levy stable noise C(i)- Increments /S.La{^t) = La{t + 
At) — La (t) of the a-stable Levy process are distributed 
according to the a-stable density with the stability index 
a. If the time step is set to At, the appropriate a-stable 
distribution is Pa(AL„; cr(Ai)i/") [11 [H, [H IH . 

The equation (|3]) is associated with the following frac- 
tional Fokker-Planck equation (FFPE) [H-Hl] 



dP{x,t) 
di 



dx 



V'{x,t) 



ga 



d\x[ 



P{x,t) (6) 



where the fractional (Riesz-Weyl) derivative is defined 

by the Fourier transform H, H^, [13, HI] J" \^^f{x) = 

— \k\°-F[f[x)]. The (space) fractional derivative in 
Eq. dni describes long jumps which are distributed ac- 
cording to the a-stable density, see Eq. ^ and Refs. [l|, 
[Sa . 1371 Wj{ . The Langevin equation Q provides a 
stochastic representation of the fractional Fokker-Planck 
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FIG. 2. Left column presents interquantile distance (go. 9 — 
qo.i, go. 8 - <7o.2, <?o.7 - go. 3 and go. 6 - go. 4, from top to bottom) 
as a function of time t. Solid lines present i^' " scaling of 
the interquantile distance, which is observed in the force free 
case. Right column demonstrates complementary cumulative 
distributions at the end of simulation. Solid lines present 
x^"' decay. Various rows correspond to different potentials: 
V{x) — (top panel), V{x) = \x\ (middle panel) and V{x) — 
X /2 (bottom panel). The stability index a and the scale 
parameter a are set to a = 0.9 and a = 1 respectively. 



equation ([6]). In other words, it describes an evolution 
of a single realization of the stochastic process {x{t)}. 
From ensemble of trajectories it is possible to analyze 
further properties of the system, in particular a proba- 
bility density P{x, t) of finding a random walker at time 
t in the neighborhood of x, see Eq. ([6|). More details on 
numerical scheme of integration of stochastic differential 
equations with respect to a-stable noises can be found in 

[ii,[il,[2i,[2i,[33. 

In the following sections properties of stationary prob- 
ability distributions for single-well systems perturbed by 
symmetric Levy noises are discussed. In the limiting 
cases of parabolic and quartic potentials the performed 
simulations corroborate earlier theoretical findings |9l.[l9l- 

El,l43. 



III. RESULTS 

Only in the limited number of special cases it is possi- 
ble to find stationary probability densities P{x) analyt- 
ically by solving Eq. ([5]). Therefore, the construction of 
stationary densities rely on numerical methods. In such 
a case, there are two frameworks possible: one option is 
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FIG. 3. Interquantile distance (go.g-go.i, <?o.8-go.2, go.7-go.3 
and go. 6 — go. 4, from top to bottom) as a function of time t for 
V{x) = IkI^"' with different values of the stabihty index a — 
{1.5, 1.1} (left column, from top to bottom) and a = {0.9, 0.8} 
(right column, from top to bottom). The scale parameter a 
is set to (T = 1. 



FIG. 4. Complementary cumulative distributions at the end 
of simulation for V{x) — \x\^'^ with different values of the 
stability index a = {1.5, 1.1} (left column, from top to bot- 
tom) and a = {0.9,0.8} (right column, from top to bottom). 
The scale parameter a is set to cr = 1. Solid lines present 



-(a-0.5) 



decay predicted by Eq. (|2ip . 



to discretize Eq. ^, see [19ll2lU42|-|45l|. alternatively one 
might use a Monte-Carlo method based on the simula- 
tion of the Langevin equation (j?]), see |25l. [26l. ISOl [33. |46| . 
Here, we rely on Monte Carlo simulations. Such a choice 
is based on statistical properties of searched densities. As 
we proceed to show, for subharmonic potentials, if sta- 
tionary states exist, the tail of the corresponding prob- 
ability density function P{x) is "fatter" than the tail of 
the corresponding a-stable density. In such a case con- 
fining the numerical solution of equation ([5]) to a finite 
interval may introduce uncontrollable errors. Therefore, 
in our simulations, we rely solely on stochastic represen- 
tation of the fractional Fokker-Planck equation (|6]) which 
is provided by Eq. ^. 

In the a = 2 case (Gaussian noise) the stationary solu- 
tion for Eq. (|6]), which in this case reduces to a "normal" 
Fokker-Planck equation, has the Boltzmann-Gibbs form 



P{x) ex exp 



V{x) 



(7) 



and exists for any power-law potential V(x) such that 



lim 



\x\- 



^Vix) 



The situation drastically changes 



for Levy noises with a < 2. Stationary states always 
exist in potentials steeper than the parabolic one |19l - 
[2l|, |4l|. However, as it will be shown this is no more the 
case for subharmonic potentials. 



A. Parabolic and quartic potentials 

There are two special cases of power-law potentials for 
which the stationary solution can be given in a closed 
form: the case of a parabolic potential (for any < a ^ 
2) and the case of a quartic potential for a = 1 and 



a = 2. These cases will serve as a benchmark for our 
further considerations. 

In order to find stationary solution of the fractional 
Fokker-Planck equation ^ we rewrite it in the Fourier 
space 



dP{k,t) 
dt 



^ U{k)P{k,t) - a°'\k\''P{k,t). 



(8) 



In Eq. ([5]) lJ{k) denotes the operator related to the 
Fourier representation of the potential V(x) |2l| . 

For the parabolic potential V{x) = x^/2 the expression 
for U{k) reads U{k) — —k-^. For symmetric a-stable 
noises, the stationary probability density fulfills [l^ l20l | 



dP{k) 
dk 



= — cr"signfc|/s| 



'Pik). 



The solution of Eq. ^ is 

P{k) = exp 



C7"\k[ 



(9) 



(10) 



i.e., the stationary solution is a symmetric Levy distribu- 
tion, see Eq. ([5]), characterized by a different value of the 
scale parameter than the noise in Eq. (j?]). The scale pa- 
rameter of the stationary density is a/a^'". For a < 2, 
the variance of the stationary solution diverges and the 
parabolic potential is not sufficient to produce bounded 
states, i.e. states characterized by a finite variance J9l.ll9l- 

mi,®. 

For the quartic potential ^(a;) = x^/4 the expression 
for U{k) reads U{k) = k-^. The stationary density 
fulfills 



d^Pjk) 
dk^ 



cr"signfc|A:|"-^P(A:). 



(11) 
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FIG. 6. Interquantile distance (go. 9 — 90.1, 90.8—90.2, 90.7 — 90.3 
and 90.6 — 90.4, from top to bottom) as a function of time t 
for V{x) = \x\ with different values of the stability index a = 
{1.7, 1.6} (left column, from top to bottom) and a = {1.5, 1.4} 
(right column, from top to bottom). The scale parameter a 
is set to cr = 1. 



FIG. 5. Quantile lines qy and |gi_y| {y = {0.9,0.8,0.7,0.6}, 
from top to bottom) as a function of 1/t (left column) and 
complementary cumulative distributions at the end of simu- 
lation (right column) for V{x) = la;|^'^ with the value of the 
stability index a = 0.7 and the scale parameter a: cr = 1 (top 
panel), a = 0.1 (middle panel) and a = 0.01 (bottom panel). 
In the right column, solid lines present x^'"^ '^' decay pre- 
dicted by Eq. ((2lJ. 



For a =: 1, the solution of Eq. dTT]) reads [11, [2C 



ai/3|^| 



y3ai/3|fc| 



i^=i(fc) = — P exp :7-^ — '- cos 

(12) 

The formula for the corresponding stationary density 
P{x) [Hill mi in the real space is 



Po.=i{x) 



1 



(13) 



TTW [1 — [x/wY + [xjwY] ' 

where w = u^/^. The stationary solution ((T^ of Eq. (ini . 
as well as the steady state (fTO]) with a < 2, is no longer 
of the Boltzmann-Gibbs type. This is typical for system 
driven by Levy white noises with the stability index a < 2 
[47| . The stationary state, sec Eq. (fT3|) . has the asymp- 
totic power-law dependence, i.e. P{x) oc \x\~^. Further- 
more, the stationary solution p^ is bimodal, with modal 
values located at x = ±a^l^l\f2 [H-iH. 



B. Subharmonic potentials 

1. Analytical results. 

Subharmonic potentials {V{x) = \x\'^ with < c < 2, 
see Fig. [T]) interpolate between the force free case (c — 0) 



and the harmonic potential (c = 2). As it was indicated 
in the previous subsection for the harmonic potential, the 
stationary state is an a-stable Levy type distribution, see 
Eqs. ([2]) and pO|) . In the force free case the stationary 
solutions do not exist. The time dependent solutions 
are just Levy stable distributions with the growing scale 
parameter cr, i.e. a{t) =at^^". This observation suggests 
that for intermediate values of the exponent c (0 < c < 2) 
there should be a transition between the situation when 
the stationary state exists and the situation when the 
stationary state is absent. 

It is possible to obtain a necessary condition for the 
existence of the stationary state using qualitative argu- 
ments based on the FFFE ^. If a stationary state for 
Eq. ([ni exists, it satisfies 



dx ^ ' d\x\°' 



P(x) 



0. 



(14) 



Assuming that tails of the distribution are given by a 
power-law, i.e. P{x) ex \x\^'^ as |x| — > 00, the first term 
in Eq. ([T4| behaves as 



Tx W-l' 



'P{x)] 



(15) 



The fractional derivative, see Eq. (|14p . is a nonlocal op- 
erator. The asymptotics of the fractional derivative of 
a probability density has a universal behavior indepen- 
dent of the particular form of this density. Indeed, its 
behavior for a; — )■ cx) is governed by the behavior of its 
Fourier-transform —\k\"P{k) for fc — >■ 0. Due to the nor- 
malization condition P{k = 0) = 1. The inverse Fourier 
transform results in 



d\x\°' 



P{x) 



-l-a 



(16) 



In order to satisfy Eq. (|14p . both terms ^TE\\ and ^TE\\ 
have to represent the same x-dependence. Consequently, 
we get c — 2 — ly = — 1 — a and 



v = c + a — I. 



(17) 



Therefore a stationary state is asymptoticahy character- 
ized by a power-law 



P{x) oc |a;| 



{c+a-l) 



(18) 



i.e. by the same expression as for the superharmonic 
potentials [2l|. The probability density function P{x) 
has to be integrable, which corresponds to i' > 1. Thus, 
the necessary condition for the existence of the steady 
state is 



c>2 



a 



a > 2 - c. 



(19) 



(20) 



For the parabolic potential (c = 2), the exponent u re- 
produces the decay of the Levy distribution {v = 1 -\~ a), 
see Eq. pUj) . For the quartic potential (c — 4) with a = 1 
we get ;/ = 4, as follows from Eq. (|13p . The complemen- 
tary cumulative density function {Fc{x) = 1 — F{x) = 
1 — /_ P{x')dx') is thus asymptotically characterized 
by a power-law decay 



Fc(x) ex x" ("+"-'). 



(21) 



Having the condition ((T^ . it is also possible to deter- 
mine the dependence of the width of the stationary dis- 
tribution w(a) on the noise strength a. Assuming that 
P{x) = w~^ f{x/w) can be expressed through a univer- 
sal function /(^) of a new, dimensionless, length vari- 
able ^ = x/w (this assumption is reasonable due to the 
scale- free nature of the power- law potential) , we change 
the length variable to ^ = x/w. For V{x) = \x\'^ (i.e. 
V'{x) (X c|a;|^~-^, x ^ 0), we get 

(22) 
The universal rescaled equation for /(^) 



|w-V(e)] + |p/(e) = o (23) 

can only be obtained if the width of distribution w(a) is 
the solution of the following equation 



w"-^ =a°'w-^-". 



(24) 



Condition ([^^ assures that prefactors of both terms in 
Eq. ((22)) scale similarly when a is changed. Thus, from 
Eq. ([M)) one obtains 



w(cr)(xa"/('=+"-2). 



(25) 



In the absence of the noise {a — >■ 0) the solution P{x) 
corresponds to the particle's position in the minimum of 
the potential, P{x) = 6{x), i.e. w{a) = 0. For c + a < 2 
the exponent in Eq. (j25p is negative, and the distribu- 
tion's width 'w{<7) diverges for ct — >■ which corresponds 
to a non-physical situation. Consequently, no station- 
ary state is possible. On the contrary, for c -I- a > 2 
the exponent in Eq. (|25p is positive and the distribution 
width is the increasing function of the scale parameter a. 
Therefore, the necessary condition of the existence of the 
steady state, see Eq. (fT9|) . is confirmed by the condition 
(HSl). 

Note that the cases of parabolic and quartic potentials 
confirm Eq. ([25]) . For the parabolic potential (c = 2) 
we get w{<7) oc a, see Eq. ([T0|) . In the case of quartic 
potential (c — 4) with a = 1 we get w{a) (x a^^^, see 
Eq. (|T3)) . Moreover, it is also possible to discuss how 
does the relaxation time to the equilibrium depends on 
the scale parameter a. To do this, we turn to Eq. ([6]) and 
rescale both the length (^ = x/w) and the time variables 
(r = i/ro(cr)), so that the rescaled solution to Eq. ([6]) is 
a universal function. The rescaling of t corresponds to 



To(cr) OC [w{a)Y 



a{3-c)/{c+a~2) 



(26) 



It indicates that the dependence on the scale parameter 
a for the stability index a close to 2 — c gets extremely 
sharp. In other words, when a approaches 2 — c (from 
above) the relaxation time grows and diverges at a = 
2-c. 

The condition (fT9|) for the existence of the steady state 
can be corroborated by the following qualitative argu- 
ment based on the decomposition of the noise into a 
"background" noise with finite variance, and "bursts" of 
arbitrary amplitude [4^, |40| . 

Let us consider the discretized Langevin scheme with 
the time step of integration Ai and the characteristic am- 
plitude of the Levy jumps e = cr(At)^/". Additionally, let 
us introduce a cutoff level A, so that only the jumps with 
the amplitude a > A are considered as bursts. The over- 
all situation is then considered as pertinent to a pertur- 
bation of the equilibrium distribution of particles in the 
potential V{x) under the action of the background noise 
(i.e. the Boltzmann distribution in the potential V{x) at 
the temperature kT ~ A) by these distinct bursts. The 
width of this distribution W depends on the cutoff level 
A and on the type of the potential. Moreover, we assume 
that only very large bursts, with a ^ A, can be respon- 
sible for the absence of a stationary distribution. There- 
fore, the finite width of the background distribution can 
be neglected when considering the action of such bursts. 
The distribution of the burst amplitudes is then given by 
the asymptotic behavior of tails of the symmetric Levy 
distribution 



p{a) 



a"' At 



x 



l + Q 



a; 



l+a 



(27) 



A burst of amplitude a brings the particle to a position 
a; « a, and the typical time necessary to return to the 



interval of the characteristic width W is given by the 
solution of the ordinary equation of motion 



x{t) = -V'{x). 



(28) 



For c 7^ 2, the return time T(a) into the equilibrium 
domain, i.e. |a;| < W, from the initial position a; = a is 
given by 



T{a) 



1 



c(2 - c) 



W^ 



(29) 



For subharmonic potentials, < c < 2, the return time 
T{a) is dominated by the initial position a. The lead- 
ing contribution to T{a) is provided by the first term of 
Eq. (gH), i.e. a'^-''/[c{2 ~ c)\. Therefore, instead of the 
return time to the equilibrium domain it is possible to 
consider the return time to the origin, which is finite and 
slightly larger than the return time to the characteristic 
domain. 

It is reasonable to assume, that the steady state does 
not exist, if the time T{a) to return from the point a 
is larger than the typical time between two bursts of the 
amplitude a or larger, and might exist in the opposite sit- 
uation. The probability to have a burst of the amplitude 
a or larger is 

P{a) = 2f p(a')da'~ (-)", (30) 

so that the typical time between two such events is 

r~Ai(-)". (31) 

Comparing T with T{a) we conclude that the inequality 
T{a) < T applies and the stationary states arc possible 
for 



c> 2 



(32) 



For superharmonic potentials, c > 2, the situation is 
inverse to the situation for subharmonic potentials. The 
larger the initial displacement a is, the faster is the return 
to the W^-domain. For the increasing initial displacement 
a, the typical return time essentially tends to a constant 
value W^~'^/[c(c— 2)], while the time between two subse- 
quent bursts grows with their amplitude a. Comparing T 
with T{a), we conclude that stationary states for super- 
harmonic potentials exist for every value of the stability 
index a. The harmonic potential, c = 2, is the limiting 
case for which the return time shows not a power-law but 
a logarithmic dependence. In such a case, the very same 
arguments proofing existence of stationary states hold. 



2. Numerical results. 

Numerical results were obtained by Monte Carlo sim- 
ulations of Eq. (m with the time step of integration 



Ai = 10 ^ and averaged over A^ = 10^ realizations 
[3 m, m [H, [33. Initially, at time i = a test parti- 
cle was located at the origin, i.e. .t(0) = 0. Due to the 
symmetry of the noise and of the potential, stationary 
states are symmetric with respect to the origin; the me- 
dian and modal values of stationary densities arc located 
at the origin. 

Our test of stationarity is based on the quantiles of dis- 
tributions. Quantiles of stationary distributions remain 
constant. Consequently, plotting quantile values as func- 
tions of time (quantile lines) gives us the key to analyze 
stationarity. More precisely, if the quantile lines are par- 
allel to the abscissa, it means that the stationary state 
has been reached. 

In order to verify performance of the test based on 
quantile lines we have constructed quantile lines for Levy 
fiights in the parabolic potential. In such a case the sta- 
tionary state exists and is given by the a-stable density, 
sec Eq. (|10p . Therefore, quantile lines should be paral- 
lel to the abscissa. Analogously, interquantile distance 
should not change over time. This behavior is indeed 
observed, e.g. in the bottom panel of Fig. [5J 

Furthermore, Fig. [5] compares interquantile distance 
(left column) and complementary cumulative distribu- 
tions Fc{x,t) = 1 - F{x,t) = 1 - J^^ P{x',t)dx' (right 
column) for the stability index a = 0.9 and various po- 
tentials: V{x) = (top panel), V{x) = \x\ (middle panel) 
and V{x) ~ x'^ /2 (bottom panel). From these cases only 
the V{x) — .T^/2 case corresponds to the situation when 
stationary state exists. Top panel of Fig. [3] corresponds 
to the force free case in which no stationary density ex- 
ists. In such a case, the interquantile distance grows like 
^i/a rpj^g middle panel of Fig. [5] corresponds to the in- 
termediate case, i.e. V{x) = \x\. Numerical simulations 
indicate that for V{x) = \x\ and a = 0.9 no stationary 
state exists as well, what is coherent with analytical ar- 
guments given above. Note that in all three cases the 
decay of the probability densities follows the one of the 
noise, albeit on different reasons. For free motion, this 
decay follows the one of the noise due to the stability of 
the process. For the case V{x) = \x\, where according to 
analytical result the stationary state is absent, numerical 
simulations suggest that the tail of the distribution func- 
tion is still dominated by the properties of the noise, and 
follows Fc{x) (X I'"" pattern. Finally, for the parabolic 
potential it is essentially given by Fc{x) oc x~'^^+"~^-', i.e. 
again by Fc{x) ex x^" for c = 2. 

Consecutive figures present results for subharmonic po- 
tentials. Figs. [SHSl present results for V{x) = |.t|^'^. 
Figs. [SHU present: interquantile distance as a function 
of time t (Fig. ^ and complementary cumulative dis- 
tributions {Fc{x) = 1 — F{x)) at the end of simulation 
(Fig. S]). Various panels present result for various val- 
ues of the stability index a: a = {1.5, 1.1} (left column, 
from top to bottom) and a = {0.9, 0.8} (right column, 
from top to bottom). Finally, Fig. [5] presents quantile 
lines (left column) and complementary cumulative dis- 
tributions (right column) for a = 0.7 with various scale 
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FIG. 7. Complementary cumulative distributions at the end 
of simulation for V{x) = \x\ with different values of the sta- 
bility index a = {1.7, 1.6} (left column, from top to bottom) 
and Q = {1.5, 1.4} (right column, from top to bottom). The 
scale parameter a is set to cr = 1. Solid lines present a;"'""^' 
decay predicted by Eq. 1)21^ . 



FIG. 8. Quantile lines qy and gi_y (y = {0.9, 0.8, 0.7, 0.6}, 
from top to bottom) as a function of 1/t (left column) and 
complementary cumulative distributions at the end of simu- 
lation (right column) for V{x) — \x\ with the value of the 
stability index a — 1.3 and the scale parameter a: a = 1 (top 
panel) and a — 0.1 (bottom panel). In the right column, solid 
lines present x~^°'^^' decay predicted by Eq. (|2ip . 



parameters a: a = 1 (top panel), cr = 0.1 (middle panel) 
and a = 0.01 (bottom panel). 

Figures [5H5] present results for V{x) = \x\. Figs. [5]- 
[7] present: interquantile distance as a function of time 
t (Fig. [6]) and complementary cumulative distributions 
iFc{x) ~ 1— F(x)) at the end of simulation (Fig. [7]). Vari- 
ous panels present result for various values of the stability 
index a: a = {1.7, 1.6} (left column, from top to bottom) 
and a = {1.5, 1.4} (right column, from top to bottom). 
Finally, Fig. [8] presents quantile lines (left column) and 
complementary cumulative distributions (right column) 
for a = 1.3 with various scale parameters cr: cr = 1 (top 
panel) and cr = 0.1 (bottom panel). 

The stability index a in Figs. [SHSl is chosen in such a 
"way that Eq. p9|) holds true. Even in this case the nu- 
merical verification "whether stationary states exist or do 
not exist, due to possible slow convergence to a station- 
ary state, is not trivial, see discussion of Figs. [5] and [5] 
below. 

The probability densities for the system described by 
Eqs. (HI and (O arc symmetric with respect to the origin. 
Therefore the quantile lines for qy and |9i-y| overlap and 
are not distinguishable in left panels of Figs. [5]and|Sl 



For V{x) 



11.5 



Fig. [3] indicates that for sufhciently 



large values of the stability index a stationary densities 
exist. Although all the values of stability indices shown 
correspond to the domain where stationary states are pre- 
dicted to exist, the convergence to these states becomes 
slower when a approaches 2 — c = 0.5 for V{x) = \x\^^^. 
Fig. [5] shows the very similar behavior for V{x) = |a;|. 
Here again the relaxation time to the equilibrium gets 
larger when a approaches 2 — c = 1. 

Figurcs|3]and[7]present complementary cumulative dis- 
tributions in the situation when Eq. (J19p is valid. There- 
fore, according to Eq. (j2ip . complementary cumulative 



distributions should asymptotically decay as a;~'^+"^^^ 
Indeed, in Figs. H] and [7] the decay predicted by Eq. (|2T|) 
is observed. 

Figures [5] (a = 0.7 and V{x) = \x\^-^) and[a(a = 1.3 
and V{x) = |a;|) present quantile lines and complemen- 
tary cumulative distributions for different values of the 
scale parameter a. The quantile lines are plotted as func- 
tions of 1/i in order to elucidate convergence. Thus, in 
Fig-Elthe convergence of 50.9 is not achieved over the time 
of simulation for cr = 1, but is approached for a ~ 0.01. 
Correspondingly, the tail of the complementary cumula- 
tive density differs from theoretical predictions for ct = 1 
and approaches its theoretical asymptotics for a = 0.01. 
The same is true for the data in Fig. [H where the con- 
vergence is achieved already for a = 0.1. This indicates 
that the problem of the slow convergence to the station- 
ary state has to be taken seriously. 



IV. SUMMARY AND CONCLUSIONS 

We considered the problem of existence of stationary 
states in power-law subharmonic potentials V{x) = \x\'^, 
c < 2, under the action of Levy-stable noise character- 
ized by the stability index a. The scaling analysis of the 
fractional Fokker-Planck equation leads us to the conclu- 
sion that such states exist if a > 2 — c. This conclusion 
is corroborated by an alternative argument based on the 
decomposition of the Levy noise. For subharmonic po- 
tentials, the probability density functions of stationary 
states are characterized by the asymptotic power-law de- 
cay P(x) oc \x\~'' = |x|~''^"''"~-^' for \x\ — )• 00, which is 
slower than decay of the tails of the corresponding Levy 
distributions. Monte Carlo simulations of the Langcvin 



equation confirm these analytical findings, at least for a 
not too close to 2 — c. The convergence to the stationary 
state for a approaching 2 — c is very slow, which poses 
considerable numerical problems, since both long times 
and very small values of scaling parameter a are required. 
However, for a small enough the tendency to converge is 
still observable. 
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